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ABSTRACT 

We test the consistency of estimates of the non- linear coupling constant /nl using non- 
Gaussian CMB maps generated by the method described in (Liguori, Matarrese and 
Moscardini 2003). This procedure to obtain non-Gaussian maps differs significantly 
from the method used in previous works on estimation of /nl- Nevertheless, using 
spherical wavelets, we find results in very good agreement with (Mukherjee and Wang 
2004), showing that the two ways of generating primordial non-Gaussian maps give 
equivalent results. Moreover, we introduce a new method for estimating the non- 
linear coupling constant from CMB observations by using the local curvature of the 
temperature fluctuation field. We present both Bayesian credible regions (assuming 
a flat prior) and proper (frequentist) confidence intervals on /nl, and discuss the 
relation between the two approaches. The Bayesian approach tends to yield lower 
error bars than the frequentist approach, suggesting that a careful analysis of the 
different interpretations is needed. Using this method, we estimate /nl = — 10^ 2 60 a ^ 
the 2a level (Bayesian) and /nl = — 10^270 (frequentist). Moreover, we find that the 
wavelet and the local curvature approaches, which provide similar error bars, yield 
approximately uncorrelated estimates of Jnl and therefore, as advocated in (Cabella 
et al. 2004), the estimates may be combined to reduce the error bars. In this way, we 
obtain Jml — —5 ± 85 and /jvx = —5 ± 175 at the 1<7 and 2a level respectively using 
the frequentist approach. 

Key words: cosmic microwave background - cosmology: theory - methods: numerical 
- methods: statistical - cosmology: observations 



1 INTRODUCTION 

Inflation is the standard paradigm for providing the initial conditions for structure formation and Cosmic Microwave Back- 
ground (CMB) anisotropy generation. In the inflationary picture, primordial adiabatic perturbations arise from quantum 
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fluctuations of the inflaton scalar field which drives the accelerated Universe expansion. In the simplest models, the infla- 
ton is assumed to have a shallow potential, thereby leading to a slow rolling of this field down its potential. The flatness 
of the potential implies that intrinsic non-linear (hence non-Gaussian) effects during slow-roll inflation are tiny, although 
non-zero and calculable l |Falk, B angarajan fc Srednicki 19 93 Gan gui et al. 1994| |Lesgourgues,"P olarski fc Starobinsky 1 997 
|Wang fc Kamio nkowski 2000 Acquav iva et al. 2 003 Maldaccna 2002). To quantitatively describe the theoretical findings in 
this framework, let us introduce a useful parameterisation of non-Gaussianity according to which the primordial gravitational 
potential $ is given by a linear Gaussian term </>g, plus a quadratic contribution, as follows (e.g. llVerde et al. 2 000)'): 

$(x) = fcj(x) + /nl4(x) (1) 

(up to a constant offset, which only affects the monopole contribution), where the dimensionless parameter /nl sets the 
strength of non-Gaussianity. The above mentioned calculation of the amount of non-Gaussianity during single-field inflation 
leads to typical values /nl ~ 10 -2 , much too low to be observable in CMB experiments. However, non-linear gravitational 
corrections after inflation unavoidably and significantly enhance the non-Gaussianity level, leading to values of /nl ~ 1, 
almost independent of the detailed inflation dynamics ( |Bartolo, Matarrese fc Riotto 2004a| |. An angular modulation of the 
quadratic term is also found (Ba rtolo, Matarrese fc Riotto 200 4b), so that /nl should be considered as a kernel in Fourier 
space, rather than a constant. The resulting effects in harmonic space might be used to search for signatures of inflationary 
non-Gaussianity in the CMB l |Liguori, Matarrese fc R iotto 2004). Nonetheless, owing to the large values of |/nl| considered 
here (~ 10 2 ) we will disregard this complication and assume /nl to be a constant parameter. Despite the simplicity of the in- 
flationary paradigm, the mechanism by which adiabatic (curvature) perturbations are generated is not yet fully established. In 
the standard scenario associated to single-field models of inflation, the observed density perturbations are due to fluctuations 
of the inflaton field, driving the accelerated expansion. An alternative to the standard scenario which has recently gained 
increasing attention is the curvaton mechanism ( |Enqvist fc Sloth 2002| IMoroi fc Takahashi 20011 IMoroi feT akahashi 2002 
|Lyth, Ungarelli fc Wands 2003| |Lyth fc Wands 2002 Ba rtolo, Matarrese fc Riotto 2004c| l , according to which the final cur- 
vature perturbations are produced from an initial isocurvature perturbation associated to the quantum fluctuations of a 
"light" scalar field other than the inflaton, the so-called "curvaton", whose energy density is negligible during inflation. 
Due to a non-adiabatic pressure perturbation arising in multi-fluid systems l)Mollerach 1990^ curvaton isocurvature pertur- 
bations are transformed into adiabatic ones, when the curvaton decays into radiation much after the end of inflation. An- 
other recently proposed mechanism for the generation of cosmological perturbations is the inhomogeneous reheating scenario 
(Dvali, Gruzinov & Zaldarriaga 2004 Zaldarriaga 2004 Kofman 2003 Matarrese & Riotto 2003). It acts during the reheating 
stage after inflation if super-horizon spatial fluctuations in the decay rate of the inflaton field are induced during inflation, caus- 
ing adiabatic perturbations in the final reheating temperature in different regions of the universe. An important feature of both 
the curvaton and inhomogeneous reheating scenarios is that, contrary to the single-field slow-roll models, they may naturally 
lead to high levels of non-Gaussianity. Large levels of non-Gaussianity are also predicted in a number of theoretical variants of 
the simplest inflationary models. First, generalised multi-field models can be constructed in which the final density perturbation 
is either strongly | |Salopek fc Bond 1990] |Salopek fc Bond 19911 IKof man 2003JI or mildly ( |Bartolo, Matarrese fc Riotto 200*21 
IBernardeau fc Uzan 200*21 iBc rnardea vTfc Uzan 20031 |Enqvist fc Va ihkoncn 2004) non-Gaussian, and generally characterised 
by a cross-correlated mixture of adiabatic and isocurvature perturbation modes (IWands et al. 2002JI . Values of |/nl| ~ 10 2 
are also predicted in the recently proposed ghost-inflation picture (. Arkani-Hamcd et al. 20041 . as well as in theories based on 
a Dirac-Born-Infeld (DBI)-type Lagrangian for the inflaton (Alishahiha, Sil verstein fc Tong 2004} . Quite recently, there has 
been a burst of interest for non-Gaussian perturbations of the type of Eq. £Q|. Different CMB datasets have been analysed, with 
a variety of statistical techniques (e.g. (Komats u et al. 20031 |Cayon et al . 2000 Sant os et al. 200 3 Mukh erjee fc Wang 2 004 
Gaztanaga & Wagg 2003)) with the aim of constraining /nl- In the last years some authors set increasingly stringent limits 
on the primordial non-Gaussianity level in the CMB fluctuations. Using a bispectrum analysis on the COBE DMR data 
(Komats u et al. 20021 found |/nl| < 1500. On the same data, jCayon et al. 2003} found |/nl| < 1100 using Spherical Mex- 
ican Hat Wavelets (SMHW) and (Sa ntos et al. 2003JI using the MAXIMA data set the limit on primordial non-Gaussianity 
to be |/nl| < 950. All these limits are at the la confidence level. The most stringent limit to date has been obtained by 
the WMAP team IjKomatsu et al. 2003^ : -58 < /nl < 134 at 95 % cl. Consistent results (an upper limit of /nl < 220 at 
a 2o confidence level) have been obtained from the WMAP data using SMHW (Mukherjee & Wang 2004J. It was shown in 
^Komatsu fc S pcrgcl 20011 that the minimum value of |/nl| which can be in principle detected using the angular bispectrum, 
is around 20 for WMAP, 5 for Planck and 3 for an ideal experiment, owing to the intrinsic limitations caused by cosmic 
variance. Alternative strategies, based on the multivariate empirical distribution function of the spherical harmonics of a 
CMB map, have also been proposed l |Hansen, Marinucci fc Vittorio 2003|IHansen et al. 2002fl . or measuring the trispectrum 
of the CMB IDe Troia et al. 20031 . 

The plan of the paper is as follows: in Section|2|we describe our method to produce the temperature pattern of the CMB 
in presence of primordial non-Gaussianity; Section [3] addresses statistical issues to constrain the non-linearity parameter /nl 
on the basis of WMAP data; finally, in Section^Jwe draw our conclusions. 
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2 MAP MAKING OF PRIMORDIAL NON-GAUSSIANITY 



The non-Gaussian CMB maps used in the following analysis have been generated by applying the numerical algorithm 
introduced by ( |Liguori, Matarrese fc M oscardini 2003), hereafter LMM. Here we summarise the various steps that define the 
whole procedure and refer the reader to that paper for further details. 

The starting point of the method is the simulation, directly in real space, of independent complex Gaussian variables 
nt m {r), with correlation functions 

<n^ mi (ri)nj 2m2 (r 2 )} = ^-^CJ , (2) 

where 5 is the Dirac delta function and 5j is Kronecker's delta. 

The linear potential multipoles &e m (r) (here $ L = 4>g) having the desired correlation properties, as expressed in terms 
of the primordial (i.e. unprocessed by the radiation transfer function) power-spectrum of the gravitational potential P$(k), 
can be obtained by convolving ni m {r) with suitable filter functions Wi{r, ri): 

*L(r)= J "dn r?n< m (ri)Wi(r,ri) . (3) 

The filters We(r,n) are defined as (see LMM) 



W e (r,n) = - / dkk 2 y/K(kJMkr)jt(kri) ; (4) 



here je are spherical Bessel function of order £. Notice that We(r, ri) can be pre-computed at the beginning of all simulations 
and then applied in Eq. @. One more advantage of this approach is that Wi{r,r\), at fixed r, are smooth functions of n 
which differ from zero only in a narrow region around n = r; as a consequence, the integral in Eq. ||3J can be estimated in a 
fast way by computing W((r, ri) in few points. 

At this point the values of the linear potential < I >L (r) can be recovered thanks to its expansion in spherical harmonics: 

$ L (r) = ^$LM^ m (f). (5) 

Cm 

The non-Gaussian contribution (modulo /nl) to the gravitational potential, $ NL (r) (here <1> NL = 4>q), is obtained directly in 
spherical coordinates by squaring $ L ; this is then harmonic-transformed by using the HEALPIX package ( Gorski ct al. 1998 ) 
to get <O0- 

Finally, the linear and non-linear contributions to the total CMB multipoles ai m = a^ m + /nlo^ are obtained by 
convolving each term with the real-space radiation transfer function Ae(r), 



A e (r) = - \dkk z A t {k)j t {kr) (6) 

(see LMM, for the formal derivation), 
aim = I drr 2 ^im(r)At(r), 

drr 2 <S>Z(r)A e (r) . (7) 

Notice that also the quantities At(r), which have been numerically estimated by using a modification of the CMBFAST code 
(Seljak & Zaldarriaga 1996) can be pre-computed and stored for all the simulations of a given model. 

In ( Kom atsu et al. 20031 |Komatsu, Spergel fc Wandelt 20 03 ) a different approach to produce non-Gaussian maps was 
adopted. Their starting point is the generation on a Fourier-space grid of a Gaussian field, which is then inverse-Fourier 
transformed and squared to get the non-Gaussian part of the gravitational potential in real space. The successive steps involve 
interpolation on a spherical grid, harmonic transforms and convolution with Ai(r), to obtain the Gaussian and non-Gaussian 
CMB multipole coefficients. 

3 TESTS OF NON-GAUSSIANITY 

In this section, we will estimate /nl from the WMAP data using two different approaches. The estimation procedures will 
be calibrated using non-Gaussian maps produced by the method described above. We have replicated the method applied in 
( Mukhcrjcc & Wang 2004 ) (hereafter MW) in order to check that the estimates using non-Gaussian maps generated with two 
different methods, give consistent results. We also introduce a new estimator of /nl based on the local curvature properties of 
the CMB fluctuation field. Finally, we observe that the estimates of /nl from the two methods are approximately uncorrelated; 
therefore we introduce a combined estimator which reduces the error bars by a factor of about y2. 
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3.1 Local curvature 

The local curvature test has been used in the flat limit to test the presence of non-Gaussianity due to cosmic strings 
lIDore et al. 20031 . and the extension to the spherical case (Hansen ct al. 2004 Cabella 2004), has been used to verify the 
asymmetries in the WMAP data. In this section we will test the power of the local curvature test to detect primordial non- 
Gaussianity. First we review the method. 
We consider a CMB temperature map T(9, <f>) normalised with its standard deviation a: 

T(e,<t,)^ T ^- T . (8) 

a 
The Hessian of this map can be written as 

jr^fT.ee T, e<t> \ ( Tg (T M - cotOT^/smO \ 

V 7> T. >H ) \ {T M - cot 9T A ) / smO (T w + l/2sin2(9T e )/sin 2 9 )' U 

where the comma denotes the ordinary derivative and the semicolon covariant derivative. In order to evaluate the derivatives 
it is convenient to go to harmonic-space: 

T{0, <f>),i = J2 ai m Y lm {6, 4>), t (10) 

and to use the recurrence relations (Varshal ovich et al. 1988JI : 

^-Y lrn (0,<t>) = imY lm (e,ct>) , (11) 

^Y tm (0,4>) = 1 y/l(l + 1) - m(m + l)Y lm+1 (0, </>)e- i '" - \ y/l(l + 1) - m(rn - l)Y lm _i(0, fie** ', 

twice to obtain the Hessian values in every point of the map, as suggested in (jSchmalzing & Gorski 19981. The presence of 
noise may render the derivatives unstable, but we have checked that the results in this paper do not change significantly when 
smoothing the map before performing the derivatives. 

The points of the renormalised map can be classified as: 

• hills where the eigenvalues of the Hessian are both positive, 

• lakes where the eigenvalues of the Hessian are both negative, 

• saddles where the eigenvalues have opposite sign. 

We calculate the proportions of hills, lakes and saddles above a certain level v in the normalised map. In this way, we obtain 
three functions of v which in the Gaussian case have a known functional form lIDore et al. 20031 . Deviations of these functions 
from this Gaussian expectation value can be used to detect non-Gaussianity. 
In order to constrain the /nl parameter our procedure is as follows: 

• we generate a set of primordial non-Gaussian maps with the method described in Section |2] with the WMAP best fit 
power spectrum and different /nl; 

• we convolve with a beam and add noise corresponding to the given experiment; 

• we apply the KpO galaxy and point source mask; 

• we degrade the maps to the resolution of nside = 256 where the derivatives are performed; 

• we count the densities of hills, lakes and saddles for each map and each level v with an extended mask to avoid the 
instabilities of the derivatives close to the boundaries of KpO (for details of the extensions, see ({Hansen et al. 2004)). In this 
way, we obtain the form of the hill, lake and saddle densities as a function of /nl; 

• we repeat the last three points for the data of the given experiment. 

• for each simulation as well as for the data, we construct a \ 2 as follows 

X 2 (/nl) = (x- (x)) T M-'(x- <x», (12) 



max 



where x = \h[yo),h(v\), ..., ft(!/ max ), l(i>o),l(vi), ..., l(Vmax)], h[y) and l(v) are the hill and lake densities and Vo,V\, ..., v n 
are the threshold values used, given by vq — — 3a and ^ max = 2a. The maximum threshold was determined in such a way 
as to obtain sufficient statistics for the lake proportions which go to zero at high thresholds. The covariance matrix M with 
elements Mij — (xiXj) is evaluated on the basis of 1000 Gaussian simulations. The distribution of the hill and lake densities 
at each threshold has been found to be close to Gaussian, justifying the above form of the x 2 ■ 

• The estimate of /nl is obtained for each simulation and for the data by minimising this \ with respect to /nl. 

_ 2 In 

• Bayesian error bars are derived by constructing the likelihood L(/nl) oc e x ' . Integrating this likelihood with respect 
to the parameter, we obtain the approximate Bayesian "credible regions" (see Section [3.4fl . 

• The frequentist error bars are derived by the histogram of the estimates /nl from simulations and defining the 1 and 2a 
levels as the limits within which 68% and 95% of the estimates fall. 



Primordial non-Gaussianity 5 

Note that we will perform the estimates (1) using only the diagonal part of the covariance matrix and (2) including the 
full covariance matrix. We will later show that this makes a huge difference when considering the error bars, and that care 
has to be taken when approximating the covariance matrix to be diagonal. 

The results of these simulations are shown in the left panel of Figure where we can see the effect for the different values 
of /nl on the hill density. We applied this procedure to the publicly available WMAP data 1 . We co-added the (foregrounds 
subtracted) maps from the three WMAP frequency channels Q, U and W following the procedure in IjBennett et al. 2003|l . In 
Figure|5]we show the \ (/nl) for the data around its minimum in the two cases, with and without the off-diagonal elements 
of the covariance matrix. When using only the diagonal parts of the covariance matrix we estimate /nl = 30 with a 68% 
and 95% credible region equal to —50 < /nl < 140 and —170 < /nl < 240, respectively, in good agreement with previously 
released estimates IKomatsu et al. 200 3 Mukhc rjee fc Wa ng 2004). We also include proper (frequentist) confidence intervals, 
which turn out to be —230 < /nl < 280 and —540 < /nl < 570, respectively. The relationship between the two approaches 
is discussed in subsection 13.41 When including the full covariance matrix, we estimate /nl = — 10 with a 68% and 95% 
approximate credible region equal to —120 < /nl < 120 and —270 < /nl < 260 respectively, while the frequentist constraints 
turn out to be —140 < /nl < 120 and —280 < /nl < 300. We observe that when not including the full covariance matrix, 
there is a huge difference in error bars between the Bayesian/frequentist approaches. Including the full matrix, this difference 
is smaller but persists. The much bigger frequentist error bars can be explained as follows: when assuming that the thresholds 
are uncorrelated, huge deviations from the expected value at many thresholds provide strong evidence against the model to 
be tested. On the other hand, taking into account the fact that the thresholds are highly correlated, the evidence becomes 
weaker and the expected model can still be consistent with the observations. Note further that Bayesian error bars tend to be 
lower than in the frequentist case in the above estimates from the WMAP data; the same phenomenon occurs in simulated 
maps. 

As a final remark, we note that equation il 1 21 can be exploited to implement a goodness-of-fit test for our specification of 
noise and foreground features. More precisely, we compared the best fit value x(/nl = 30) obtained for the WMAP data 
with minx(/NL) from 200 Monte Carlo simulations of non-Gaussian CMB maps with /nl = 30. When using the diagonal 
covariance matrix, the observed value corresponds to the 28% quantile, thereby suggesting good agreement between WMAP 
data and our simulated models. When we adopt the full covariance matrix, the observed value is smaller than 92% of the 
simulations, that is, the local curvature statistics on WMAP data are closer to their expected values than the great majority 
of simulated maps. Note however that we are performing four similar tests on the data (local curvature, wavelets with and 
without covariance matrix) and therefore a 92% deviation from simulations is not significant. 



3.2 Spherical Wavelets 

Wavelets are a very flexible tool used in connection with CMB data for denoising fSa nz et al. 1999"|l . extracting point sources 
jCayon et al. 2 000 Tcn orio et al. 1999ll . and detecting non-Gaussianity HUabella et al. 2004 Martinez- Gonzal ez et al. 2002 
IBarreiro fc Ho bson 200l] |Forni fc Aghanim 200l] |Starck, Aghanim ~fc Forni 2 004). To test the consistency of our maps of pri- 
mordial non-Gaussian models with those of (IKomatsu et al. 200311 we have repeated the method for estimating /nl described 
in MW. The basics steps for the test are the following (for more details, see MW): 

• we generate a set of primordial non-Gaussian maps with the method described in Section [2] with the WMAP best fit 
power spectrum and different /nl; 

• we add noise and convolve with a beam corresponding to the given experiment; 

• we apply the KpO galactic cut, leaving the point sources unmasked; 

• we degrade the maps to the resolution of nside = 256; 

• we obtain the wavelet coefficients by convolving each map with the Spherical Mexican Hat Wavelet (SMHW) given by: 



9(y,R) = 



2tyN(R) 



(I) 1 



2-1-1 e 



V\ 2 -yl/2Rl 



(13) 



using the scales R given in MW (14, 25, 50, 75, 100, 150, 200, 250, 300, 400, 500, 600, 750, 900, 1050 arcmin) ; 

• we use only the coefficients outside an extended KpO mask (for details of the extensions, see llVielva et al. 20041 and 
MW) to obtain the skewness as a function of /nl- The result is shown in the right panel of Figure 

• we finally repeat the last three steps for the data of the given experiment. 

• for each simulation as well as for the data, we construct a \ 2 as follows 

X 2 (/nl) = (x- <x)) T M- 1 (x- (x)), (14) 

1 obtainable from the LAMBDA website: http://lambda.gsfc.nasa.gov/ 
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Figure 1. On the left, density of hills as a function of the threshold v, for different values of /nl- All values are divided by the Gaussian 
expectation value. On the right, skewncss of wavelet coefficients for different scales R. The solid line represents Gaussian simulations, 
the dotted (dashed) lines represent non-Gaussian simulations for different negative (positive) values of /nl; the grey band is the la 
confidence level taking into account the noise and beam of the co-added Q+V+W map with an extended KpO mask. The values of |/nlI 
shown are 50, fOO, 200, 300, 400 and 500, showing an increasing distance from the Gaussian expectation value. The bold crosses show 
the results on the data. 
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Figure 2. The plot shows the x as a function of /nl around its minimum, using the curvature test (left) and the wavelets (right) on 
the co-added Q+W+V WMAP map. Results including the full covariancc matrix are plotted as a dotted line, whereas the solid line 
corresponds to the diagonal matrix approximation. 



where the elements of the vector x are Xi — S(Ri) where S(Ri) is the skewness of the wavelet coefficients on scale Ri. The 
covariance matrix M with elements My = (xiXj) is evaluated on the basis of 1000 Gaussian simulations. The distribution of 
the skewness at each scale is close to Gaussian, justifying the above form of the \ 2 ■ 

• For each simulations and for the data we estimate /nl by minimising this % 2 with respect to the parameter. 



Bayesian error bars are found by constructing the likelihood L(/nl) oc e 



-x'/a 



Integrating this likelihood with respect to 



the parameter, we obtain the approximate Bayesian credible regions (see Section f3. 41 . 
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• The frequentist error bars are derived by the histogram of the estimates /nl from simulations and defining the 1 and 2a 
levels as the limits within which 68% and 95% of the estimates fall. 

Again, the test was performed on the co-added (foreground cleaned) Q+V+W WMAP map. The plot of the \ 2 f° r the 
data is shown in Figure [5] We estimate the Bayesian credible region, considering only the diagonal part of the covariance 
matrix (the results presented in MW were obtained using the same approximation), for /nl to be 40 ±90 at the la level and 
40 ± 180 at the 2a level, in agreement with the values derived in MW (see again the discussion in Section|23J. Considering the 
full covariance matrix, we estimate /nl to be ± 120 at the la level and ± 240 at the 2a level As for the curvature, we have 
also calculated the frequentist confidence intervals which in the first case turn out to be quite larger, — 100 < /nl < 170 and 
—230 < /nl < 320 respectively, while considering the full covariance matrix we find for the error bars on /nl —130 < /nl < 120 
and —260 < /nl < 230. Again the Bayesian credible regions are generally smaller than the frequentist confidence levels, but 
considering the full covariance matrix, the difference is reduced. Note further that the Bayesian credible regions are slightly 
underestimated when approximating the covariance matrix to be diagonal. 

As before, we also implemented a goodness-of-fit test based on Eq. Ijl4^ : we found that the values on the WMAP 
data corresponds to the 32% and 54% quantile obtained from 200 simulations with /nl = 40 when using the diagonal and 
full covariance matrices, respectively. Again, this suggests that our specification of noise and foreground masks provides a 
reasonable approximation to the experimental settings of WMAP. 



3.3 The combined test 

In this section we apply the combined procedure introduced in fCabel la et al. 200411 to improve the constraints on /nl- First 
of all, we estimated the correlation between the two estimators /nl" and / N lT by 200 Monte Carlo simulations and found 

. - cur <■ wav, . - cur,." wav , 

(/nl /nl ) - (/nl )(/nl ) ^ 0g6 , 15 s 

-y cur -y wav ' ^ ' 

suggesting that the correlation between the estimators is less than 10%. Then, we can lower the error bars on /nl by combining 
the two statistics. More precisely, we use the estimates obtained with the full covariance matrix and we evaluate the frequentist 
confidence intervals on the following statistic: 

rwav _[_ fcur 
pcomb JNL + iNL ffc\ 

/nl = 2 ' ( 16 ) 

Using the combined test on the WMAP data, we estimate /nl = —5, with the constrains at la and 2a levels of —90 < /nl < 80 
and — 180 < /nl < 170, respectively. 



3.4 A remark on confidence intervals and credible regions on /nl 

In this subsection we present a brief discussion on the evaluation of (frequentist) confidence intervals and (Bayesian) credible 
regions on /nl- Let us denote by / NL our estimate of the non-linearity parameter obtained by minimising Eqs. 112H and II H : 
also, let us denote by £(/nl;/nl) the likelihood function. As well known, an (1 — a)-level confidence interval for /nl based 
upon the observation / NL is the set of all values /nl such that 

i(/NL| /NL)d/NL > CX.. 

I /NL — /NL I > I /j^ L — /NL I 

Note that the integral is taken with respect to the estimate /nl that is, we are viewing L(/nl; /nl) as the probability density 
of our estimator. In words, we include in the confidence interval all the values /nl such that the probability to get an estimated 
parameter as / NL or further away is at least as large as a. More clearly, a value is included provided it does not entail that 
observing what we observed is less probable than a. Of course, in the special case where the distribution of /nl is Gaussian 
with mean /nl and variance a 2 which does not depend on /nl (or, in general, is symmetric under the exchange of /nl and 
/nl), we have 

L(/nl; /nl)(^/nl (17) 

/nl - /nlI>I/^ l — /nlI 

— exp|--(/ N L -/nl) 2 |c(/nl (18) 

{-^(/NL-/NL) 2 }d/NL (19) 



/nl-/nlI>I/^ l -/nl' v2iva 



/nl-/nlI>I/^ l -/nlI v2ttct 



1 

exp 
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Figure 3. The distribution of estimates /jvl as a function of the model parameter /jvz,- The isocontours of the distribution is plotted 
in terms of percentages of the number of maps used in the simulations. The distribution is clearly skewed with respect to the diagonal 
shown for clarity. The estimates of /"Mr were obtained using the curvature method (left plot) and wavelets (right plot). 



by the symmetry of the previous expression with respect to (/nl , ,/nl ) • This justifies the common practice to derive confidence 
intervals by integrating the likelihood. Rigorously speaking, however, this is no longer justified if the integrand is not symmetric 
with respect to an exchange of /nl with /nl (for instance if a is not constant with respect to /nl). One may then try to justify 
the integration of the likelihood by a Bayesian viewpoint, by assuming a flat prior and viewing L(/nl;/nl) as a posterior 
density function. The resulting set, however, should not be labelled as confidence interval (which is a frequentist concept): it 
is a Bayesian credible region, which will depend in general on the choice of the prior. 

It is occasionally stated that this dependence is overcome by the choice of a flat prior. The latter is claimed to be non- 
informative by definition: indeed, no physicist would consider a priori equally likely that /nl lies in [—10, 10] rather than it 
is in [317,337] (say); thus a flat prior, although unphysical, is usually justified as a panacea to get objective results. This 
argument is to some extent misleading, though, as it can be shown by standard counterexamples. Take for instance Eq. Q, 
and assume for brevity's sake that /nl > (otherwise duplicate our argument by symmetry). Then Eq. Q can be rewritten 



r = (£ G (x) + (ffNL<Mx)) 2 



(20) 



that is /nl = <?nl- Of course, from the physical point of view there is no reason to prefer the alternative specification in 
Eq. Q to 1201 1 . Now let us assume we impose a "non-informative" (flat) prior on <?nl :, the posterior probability becomes 



-^expl^/NL-SNL) 2 }. 



27TCT 



exp 



{~2^ NL ~~ /nl) 2 j 



dfv 



and the credible sets are thus obviously affected, although we are working with exactly the same model as before and we 
are claiming to have used no a priori information. In short, flat priors are simply shifting the choice from the form of the 
prior probability to the form of the statistical parametrisation; the latter, moreover, is not due to physical considerations, but 
simply to computational convenience. 

With this in mind, we will now use /nl estimated from 200 simulations (with the experimental settings used above) 
to investigate whether the distribution of /nl as a function of the model /nl is symmetric with respect to an exchange of 
/nl and /nl- In Figures ©, @ we show a contour plot of this function for the curvature and wavelet test, respectively. 
When not considering the full covariance matrix, this function is markedly skewed for the wavelet test, but even more so for 
the curvature test (see fig. EJ ■ This may explain the big difference between the confidence intervals and Bayesian credible 
regions for the wavelets and the even bigger difference for the curvature test when correlations are neglected. On the other 
hand, when correlations are properly taken into account, the model fit is better and the function becomes considerably more 
symmetric (see fig. 0J. This is also reflected in the fact that the difference between Bayesian and frequentist error bars is 
markedly reduced in this case . 
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Figure 4. As Fig. [3] but by using the full correlation matrix. The skewness of the distribution is clearly reduced with respect to the 
diagonal case (see Fig. |3). The estimates of /jvjL were obtained using the curvature method (left plot) and wavelets (right plot). 



4 CONCLUSIONS 

In this paper we have used simulated CMB maps of primordial non-Gaussian models, generated with the algorithm described 
in ( |Liguori, Matarrese fc Moscardini 2003[ l to estimate the non-linear coupling parameter /nl from the WMAP data. As all 
other estimates of /nl in the literature have been based on non-Gaussian CMB maps generated by a different approach 
(Komats u et al. 2003L one of the aims of this paper was to check whether a different way of generating CMB maps with the 
same kind of primordial non-Gaussianity gives a consistent estimate of /nl. In order to perform a direct test of consistency, 
we applied the estimator of /nl using SMHW presented in | |Mukherjee fc Wa ng 2004 ) on the WMAP data. When neglecting 
the scale-scale correlations, we find that our estimate of /nl is in full agreement with theirs and thus that the two ways of 
generating primordial non-Gaussian maps give fully consistent results. We have also presented a new method to estimate 
/nl based on the local curvature of the CMB fluctuation field. Our estimate of /nl with this method on the WMAP data is 
consistent with estimates of /nl using other approaches. 

Moreover, we point out the importance of including the full covariance matrix in the \ 2 when estimating /nl with these 
two methods. We also discuss the difference between frequentist confidence intervals and Bayesian credible regions on the 
estimate of /nl • We show that the difference between these two methods for finding error bars can be huge when not including 
the full covariance matrix, particularly for the curvature estimator where correlations between thresholds are important. For 
the WMAP data, the Bayesian credible region expands whereas the frequentist confidence intervals shrink when taking into 
account correlations in the covariance matrix. We find further that the Bayesian credible regions are in general smaller than 
the frequentist confidence intervals. We conclude that care must be taken when approximating the confidence intervals on 
/nl using the integral of the likelihood with respect to the parameter. 



Including the full covariance matrix, we find with the curvature method /nl 



-lOZ 



at the 2a level using frequentist 



"10-260 us i n g Bayesian credible regions. Similarly for the wavelet method, /nl 



Oi^ with 



confidence intervals, and /nl = 

frequentist confidence intervals and /nl = O^lo- Finally, we show that the two methods provide approximately uncorrelated 
estimates of /nl; this observation naturally suggests an improved estimator which combines the two methods. With this 
combined test, we find /jvl = —5 ± 85 and }nl — — 5 ± 175 at the la and 2a level respectively. 
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